function [Y,Xf,Af] = NNFunction (X,Xi,~)
%MYNEURALNETWORKFUNCTION neural network simulation function.
%
% 
% [Y,Xf,Af] = myNeuralNetworkFunction(X,Xi,~) takes these arguments:
% 
%   X = 2xTS cell, 2 inputs over TS timsteps
%   Each X{1,ts} = 1xQ matrix, input #1 at timestep ts.
%   Each X{2,ts} = 1xQ matrix, input #2 at timestep ts.
% 
%   Xi = 2x10 cell 2, initial 10 input delay states.
%   Each Xi{1,ts} = 1xQ matrix, initial states for input #1.
%   Each Xi{2,ts} = 1xQ matrix, initial states for input #2.
% 
%   Ai = 2x0 cell 2, initial 10 layer delay states.
%   Each Ai{1,ts} = 8xQ matrix, initial states for layer #1.
%   Each Ai{2,ts} = 1xQ matrix, initial states for layer #2.
% 
% and returns:
%   Y = 1xTS cell of 2 outputs over TS timesteps.
%   Each Y{1,ts} = 1xQ matrix, output #1 at timestep ts.
% 
%   Xf = 2x10 cell 2, final 10 input delay states.
%   Each Xf{1,ts} = 1xQ matrix, final states for input #1.
%   Each Xf{2,ts} = 1xQ matrix, final states for input #2.
% 
%   Af = 2x0 cell 2, final 0 layer delay states.
%   Each Af{1ts} = 8xQ matrix, final states for layer #1.
%   Each Af{2ts} = 1xQ matrix, final states for layer #2.
% 
% where Q is number of samples (or series) and TS is the number of timesteps.

%#ok<*RPMT0>

  % ===== NEURAL NETWORK CONSTANTS =====
  
  % Input 1
  x1_step1_xoffset = 49.95;
  x1_step1_gain = 0.00293556436224864;
  x1_step1_ymin = -1;
  
  % Input 2
  x2_step1_xoffset = 49.95;
  x2_step1_gain = 0.0029247890495898;
  x2_step1_ymin = -1;
  
  % Layer 1
  b1 = [1.7753226421836283411;1.2060190066304077305;0.55700134104011167491;0.17232016405962113925;-0.70669456077451042741;1.7059182102417902449;-1.0208328105622099802;1.5335325293031143445];
  IW1_1 = [0.16093635434873101997 0.50070102258768656167 1.0117218870060902258 0.056285876932266387662 0.41455553563281455354 0.29125783071271293956 0.40300657337804624003 -0.48278382062523300666 -0.50021490324304307418 -0.046010396828445696005;-0.53623872498648950113 0.29089087814605341098 -0.41731448239142454604 0.28816021535320385238 -0.78768120053554380178 -0.28067000734868124967 0.89280448204958717806 -0.70844779986838313235 0.66808314182243111468 -0.17756604892417363062;-0.61747141909604918997 0.072032489359829518505 -0.021637111758942017409 -0.20315425109473400633 -0.2874324858449580633 -0.29976258956718115911 -0.074223579321382002205 0.10548969362192239796 0.58645467533626216561 -0.50652206384249054416;0.090051102331473384388 0.20702275242301601543 -0.17044546959951131782 -0.059650292863389846321 0.59707361934398528458 -0.014671901665977937074 -0.097900898499317323687 0.27589399819657139945 -0.45812705994100932916 1.1469352441992173297;-0.78097494635131758933 -0.68336990184087842426 1.0440012414771611127 0.42999240347106754134 -0.016733175680228137144 0.011540974297546583577 -0.14084202953828359184 0.45530135804025761992 0.30488957562427315695 0.23685774014716742819;0.9638802993668182939 -0.11973832226137651191 -0.34511312141095651373 -0.03086674960184669983 0.10933803912270442793 1.3071661754589807281 0.83579258842047943201 -0.44665026304820010239 -0.51885961121422030651 0.008778374585347310799;-0.46659639230384664099 0.98622573343648634037 0.057457698315120292465 0.17178069663469802686 -1.3056552755365649077 0.52878757472086101998 -0.15324124232604094398 -1.005065241964211431 0.76349145826684672311 1.0586438035029097637;0.73652679757427585994 0.60308000156740793063 -0.61598202216050945701 -0.36888561744286313404 -0.23806161175493978366 0.51921162792003794362 0.4000099364978637051 -0.17329027959130682834 -0.026297785831993863298 1.1951561796627478973];
  IW1_2 = [-0.83641980512160507555 -0.073230490964933198095 -0.1883088249612873688 0.040496453470726258284 0.50165094279284072609 0.61238351041442107281 0.56108630923448932215 0.40237969229472569932 0.25609734139171985712 -0.26961155865108721308;-0.81730679753840296264 0.033837355154797250922 0.39905464738693324778 -0.52139723418553196588 0.33152705966647727509 -0.32066118797552867203 -0.72561120404417334129 1.0381080410381917645 -0.9340145570006584963 -0.19500381071308087177;-0.22393267083752094182 0.13560335177230772863 -0.097848602210608995322 0.27730977984631099975 0.41331960273480838897 -0.12025253722767043074 -0.21576068118308686095 0.2501880215846626343 -0.22255648969969485385 0.075083778871587736337;-1.8402416169322071138 0.77250358159968046667 0.006093525577395190497 -0.26203108824186871484 -0.21543170973757236331 -0.21478870313284120774 0.56604190050702396508 0.17132983870006901483 -0.22757491473788934999 -0.57605192538333982455;-0.7938696866564186605 -0.50329955636330347968 -0.031415923185297053855 0.28364351398692794204 0.89774599390993803372 -0.47953186504996808681 -0.58649408343059650228 0.44523735531080277372 -0.13199763201946090874 0.96222404289702434532;-1.9649946110141331967 1.0604888102164762298 -0.39073216977506597392 -0.57457882692695849691 -0.17535986226242880615 0.50790880012482164663 1.1073898376153508138 0.16094787973418206439 -1.3978978943758146869 -0.34397940283626438074;-3.5383501679607363144 -0.54127109745439783151 0.83461010414688885373 0.35124575694442067286 0.72881136584846040183 -1.0841436561135338046 0.44274482546328131072 0.64665132352265486393 -0.87850594456356201345 0.5350811508689523599;-0.3136779627818777616 0.44053722004317136252 -0.19008588388298808258 0.11514426660885822284 0.0059439029985793329774 -0.063935123636786356238 0.22962863293407376197 0.31288635002697512411 -0.58118875687768101201 -0.81514642981249885167];
  
  % Layer 2
  b2 = 0.52860785049250891277;
  LW2_1 = [0.030943129571076838369 0.14055192463124893676 -1.0706747231671522869 -0.49217815898028044463 0.073317438421473327126 -0.58707593704201865847 -0.022344814478422243831 0.57786211274042553754];
  
  % Output 1
  y1_step1_ymin = -1;
  y1_step1_gain = 0.0029247890495898;
  y1_step1_xoffset = 49.95;
  
  % ===== SIMULATION ========
  
  % Format Input Arguments
  isCellX = iscell(X);
  if ~isCellX, X = {X}; end;
  if (nargin < 2), error('Initial input states Xi argument needed.'); end
  
  % Dimensions
  TS = size(X,2); % timesteps
  if ~isempty(X)
    Q = size(X{1},2); % samples/series
  elseif ~isempty(Xi)
    Q = size(Xi{1},2);
  else
    Q = 0;
  end
  
  % Input 1 Delay States
  Xd1 = cell(1,11);
  for ts=1:10
    Xd1{ts} = mapminmax_apply(Xi{1,ts},x1_step1_gain,x1_step1_xoffset,x1_step1_ymin);
  end
  
  % Input 2 Delay States
  Xd2 = cell(1,11);
  for ts=1:10
    Xd2{ts} = mapminmax_apply(Xi{2,ts},x2_step1_gain,x2_step1_xoffset,x2_step1_ymin);
  end
  
  % Allocate Outputs
  Y = cell(1,TS);
  
  % Time loop
  for ts=1:TS
  
    % Rotating delay state position
    xdts = mod(ts+9,11)+1;
  
    % Input 1
    Xd1{xdts} = mapminmax_apply(X{1,ts},x1_step1_gain,x1_step1_xoffset,x1_step1_ymin);
    
    % Input 2
    Xd2{xdts} = mapminmax_apply(X{2,ts},x2_step1_gain,x2_step1_xoffset,x2_step1_ymin);
    
    % Layer 1
    tapdelay1 = cat(1,Xd1{mod(xdts-[1 2 3 4 5 6 7 8 9 10]-1,11)+1});
    tapdelay2 = cat(1,Xd2{mod(xdts-[1 2 3 4 5 6 7 8 9 10]-1,11)+1});
    a1 = tansig_apply(repmat(b1,1,Q) + IW1_1*tapdelay1 + IW1_2*tapdelay2);
    
    % Layer 2
    a2 = repmat(b2,1,Q) + LW2_1*a1;
    
    % Output 1
    Y{1,ts} = mapminmax_reverse(a2,y1_step1_gain,y1_step1_xoffset,y1_step1_ymin);
  end
  
  % Final Delay States
  finalxts = TS+(1: 10);
  xits = finalxts(finalxts<=10);
  xts = finalxts(finalxts>10)-10;
  Xf = [Xi(:,xits) X(:,xts)];
  Af = cell(2,0);
  
  % Format Output Arguments
  if ~isCellX, Y = cell2mat(Y); end
end

% ===== MODULE FUNCTIONS ========

% Map Minimum and Maximum Input Processing Function
function y = mapminmax_apply(x,settings_gain,settings_xoffset,settings_ymin)
  y = bsxfun(@minus,x,settings_xoffset);
  y = bsxfun(@times,y,settings_gain);
  y = bsxfun(@plus,y,settings_ymin);
end

% Sigmoid Symmetric Transfer Function
function a = tansig_apply(n)
  a = 2 ./ (1 + exp(-2*n)) - 1;
end

% Map Minimum and Maximum Output Reverse-Processing Function
function x = mapminmax_reverse(y,settings_gain,settings_xoffset,settings_ymin)
  x = bsxfun(@minus,y,settings_ymin);
  x = bsxfun(@rdivide,x,settings_gain);
  x = bsxfun(@plus,x,settings_xoffset);
end
